###################################################################################################

#Figure3_PanelA
#This figure plots the evolution of the average value of commodities disbursed from January to November of 2017 by treatment and control
###################################################################################################

### Setup #################################################################################

#DEFAULT COLORS
fColor <- "black"
mColor <- "gray60"
lightColor <- "#97B8DE"
darkColor <- "#4669BD"



ration.value_total <- c("totalr","totalnr")
ration.name_total <- c("total reconciled","total non-reconciled")


month.list <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun","Jul", "Aug","Sep","Oct", "Nov")

rect_r <- data.frame(xmin=7, xmax=10, 
                     ymin=-Inf, ymax=Inf)

fontSize <- 24
theme_set(theme_gray(base_size = fontSize))


ration.value_total <- c("totalr","totalnr")
ration.name_total  <-c("total reconciled","total non-reconciled")




##############################################################
#Reconciled
#############################################################


tmp.outcome <- "value"

r <- 1

tmp.ration <- ration.value_total[r]
tmp.ration1 <- ration.name_total[r]

tmp.data.c <- foreign::read.dta(paste(AdminDataDir,tmp.ration,"_",tmp.outcome,"_meandisbursement_reconciliation_time_series_treatment0.dta",sep = ""))
tmp.data.t <- foreign::read.dta(paste(AdminDataDir,tmp.ration,"_",tmp.outcome,"_meandisbursement_reconciliation_time_series_treatment1.dta",sep = ""))



#reshape the dataset to include values for month and whether predicted/observed
meltdf.c <- reshape2::melt(tmp.data.c[,c("month",paste("dis_",tmp.outcome,"_perRC_",tmp.ration,sep = ""),"yhat")],id="month")
meltdf.c$group <- NA
meltdf.c$group[meltdf.c$variable == paste("dis_",tmp.outcome,"_perRC_",tmp.ration,sep = "")] <- "Observed_control"
meltdf.c$group[meltdf.c$variable =="yhat"] <- "Fitted value_control"


meltdf.t <- reshape2::melt(tmp.data.t[,c("month",paste("dis_",tmp.outcome,"_perRC_",tmp.ration,sep = ""),"yhat")],id="month")
meltdf.t$group = NA
meltdf.t$group[meltdf.t$variable == paste("dis_",tmp.outcome,"_perRC_",tmp.ration,sep = "")] <- "Observed_treatment"
meltdf.t$group[meltdf.t$variable =="yhat"] <- "Fitted value_treatment"

finaldf <- rbind( meltdf.c,  meltdf.t)


#reshape to get se
meltdf.se.t <- reshape2::melt(tmp.data.t[,c("month",paste("dis_",tmp.outcome,"_perRC_",tmp.ration,"_se",sep = ""),"yhat_se")],id="month")
meltdf.se.t$se <- 0
meltdf.se.t$se[meltdf.se.t$variable =="yhat_se"] = meltdf.se.t$value[meltdf.se.t$variable =="yhat_se"]  #leave the rest zero, we don't need to plot se of observed
meltdf.se.t$value <- NULL
meltdf.se.t$group <- NA
meltdf.se.t$group[meltdf.se.t$variable == paste("dis_",tmp.outcome,"_perRC_",tmp.ration,"_se",sep = "")] <- "Observed_treatment"
meltdf.se.t$group[meltdf.se.t$variable =="yhat_se"] <-"Fitted value_treatment"
meltdf.se.t$variable <- NULL

meltdf.se.c <- reshape2::melt(tmp.data.c[,c("month",paste("dis_",tmp.outcome,"_perRC_",tmp.ration,"_se",sep = ""),"yhat_se")],id="month")
meltdf.se.c$se <- 0
meltdf.se.c$se[meltdf.se.c$variable =="yhat_se"] = meltdf.se.c$value[meltdf.se.c$variable =="yhat_se"]  #leave the rest zero, we don't need to plot se of observed
meltdf.se.c$value <- NULL
meltdf.se.c$group <- NA
meltdf.se.c$group[meltdf.se.c$variable == paste("dis_",tmp.outcome,"_perRC_",tmp.ration,"_se",sep = "")] <- "Observed_control"
meltdf.se.c$group[meltdf.se.c$variable =="yhat_se"] <- "Fitted value_control"
meltdf.se.c$variable <- NULL


finaldf.se <- rbind(meltdf.se.c,  meltdf.se.t)

plot.dat1 <- merge(finaldf,finaldf.se, by = c("month","group"))
plot.dat1$lwr <- plot.dat1$value - 1.96*plot.dat1$se
plot.dat1$upr <- plot.dat1$value + 1.96*plot.dat1$se

Upper <- max(plot.dat1$upr)


##############################################################
#Unreconciled
#############################################################
r <- 2

tmp.ration <- ration.value_total[r]
tmp.ration1 <- ration.name_total[r]

tmp.data.c <- foreign::read.dta(paste(AdminDataDir,tmp.ration,"_",tmp.outcome,"_meandisbursement_reconciliation_time_series_treatment0.dta",sep = ""))
tmp.data.t <- foreign::read.dta(paste(AdminDataDir,tmp.ration,"_",tmp.outcome,"_meandisbursement_reconciliation_time_series_treatment1.dta",sep = ""))



#reshape the dataset to include values for month and whether predicted/observed
meltdf.c <- reshape2::melt(tmp.data.c[,c("month",paste("dis_",tmp.outcome,"_perRC_",tmp.ration,sep = ""),"yhat")],id="month")
meltdf.c$group <- NA
meltdf.c$group[meltdf.c$variable == paste("dis_",tmp.outcome,"_perRC_",tmp.ration,sep = "")] <- "Observed_control"
meltdf.c$group[meltdf.c$variable =="yhat"] <- "Fitted value_control"


meltdf.t <- reshape2::melt(tmp.data.t[,c("month",paste("dis_",tmp.outcome,"_perRC_",tmp.ration,sep = ""),"yhat")],id="month")
meltdf.t$group <- NA
meltdf.t$group[meltdf.t$variable == paste("dis_",tmp.outcome,"_perRC_",tmp.ration,sep = "")] <- "Observed_treatment"
meltdf.t$group[meltdf.t$variable =="yhat"] <- "Fitted value_treatment"

finaldf <- rbind( meltdf.c,  meltdf.t)


#reshape to get se
meltdf.se.t <- reshape2::melt(tmp.data.t[,c("month",paste("dis_",tmp.outcome,"_perRC_",tmp.ration,"_se",sep = ""),"yhat_se")],id="month")
meltdf.se.t$se <- 0
meltdf.se.t$se[meltdf.se.t$variable =="yhat_se"] <- meltdf.se.t$value[meltdf.se.t$variable =="yhat_se"]  #leave the rest zero, we don't need to plot se of observed
meltdf.se.t$value <- NULL
meltdf.se.t$group <- NA
meltdf.se.t$group[meltdf.se.t$variable == paste("dis_",tmp.outcome,"_perRC_",tmp.ration,"_se",sep = "")] <- "Observed_treatment"
meltdf.se.t$group[meltdf.se.t$variable =="yhat_se"] <- "Fitted value_treatment"
meltdf.se.t$variable <- NULL

meltdf.se.c <- reshape2::melt(tmp.data.c[,c("month",paste("dis_",tmp.outcome,"_perRC_",tmp.ration,"_se",sep = ""),"yhat_se")],id="month")
meltdf.se.c$se <- 0
meltdf.se.c$se[meltdf.se.c$variable =="yhat_se"] <- meltdf.se.c$value[meltdf.se.c$variable =="yhat_se"]  #leave the rest zero, we don't need to plot se of observed
meltdf.se.c$value <- NULL
meltdf.se.c$group <- NA
meltdf.se.c$group[meltdf.se.c$variable == paste("dis_",tmp.outcome,"_perRC_",tmp.ration,"_se",sep = "")] <- "Observed_control"
meltdf.se.c$group[meltdf.se.c$variable =="yhat_se"] <- "Fitted value_control"
meltdf.se.c$variable <- NULL


finaldf.se<-rbind(meltdf.se.c,  meltdf.se.t)

plot.dat2 <- merge(finaldf,finaldf.se, by = c("month","group"))
plot.dat2$lwr <- plot.dat2$value - 1.96*plot.dat2$se
plot.dat2$upr <- plot.dat2$value + 1.96*plot.dat2$se

Lower <- min(plot.dat2$lwr)



##############################################################
# Combine
#############################################################
plot.dat2$group <- paste(plot.dat2$group,"non-reconciled")
plot.dat1$group <- paste(plot.dat1$group,"reconciled")
plot.dat <- rbind(plot.dat1,plot.dat2)





##############################################################
# Plot
#############################################################
ggplot(plot.dat, aes(x=month, y=value, group = group, colour= group, linetype =group)) +
  geom_line(size= 1.5)+ 
  geom_ribbon(aes(ymin=lwr,ymax=upr), alpha = 0.1, colour=NA) + 
  scale_colour_manual(values=c(fColor,darkColor,fColor,darkColor,mColor,lightColor,mColor,lightColor)) +
  scale_linetype_manual(values=c("dashed","dashed","solid","solid", "dotted","dotted","dotted","dotted"))+
  xlab("Month")+  
  ylab("Average disbursement value per rationcard")+ 
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16),
        legend.position="none",
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  scale_y_continuous(limit=c(Lower  ,550), breaks = seq(0,550,100)) +
  scale_x_continuous(limit=c(1,11), breaks = seq(1,11,1), labels = month.list) +
  geom_rect(data=rect_r, aes(xmin=xmin, xmax=xmax, 
                             ymin=ymin, ymax=ymax), 
            color="gray100", alpha=0.05, inherit.aes = F,
            linetype = 0)



ggsave(file=paste(OutputDir, "Figure3_PanelA.pdf", sep=""), width=12, height=6)







